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Abstract 

The Broad Histogram is a method allowing the direct calculation of the 
energy degeneracy g(E). This quantity is independent of thermodynamic 
concepts such as thermal equilibrium. It only depends on the distribution of 
allowed (micro) states along the energy axis, but not on the energy changes 
between the system and its environment. Once one has obtained g(E), no 
further effort is needed in order to consider different environment conditions, 
for instance, different temperatures, for the same system. 

The method is based on the exact relation between g(E) and the mi- 
crocanonical averages of certain macroscopic quantities N up and N dn . For 
an application to a particular problem, one needs to choose an adequate in- 
strument in order to determine the averages < N up (E) > and < N dn (E) >, 
as functions of energy. Replacing the usual fixed-temperature canonical by 
the fixed-energy micro canonical ensemble, new subtle concepts emerge. The 
temperature, for instance, is no longer an external parameter controlled by 
the user, all canonical averages being functions of this parameter. Instead, 
the microcanonical temperature T m (E) is a function of energy defined from 
g(E) itself, being thus an internal (environment independent) characteristic 
of the system. Accordingly, all microcanonical averages are functions of E. 

The present text is an overview of the method. Some features of the 
microcanonical ensemble are also discussed, as well as some clues towards the 
definition of efficient Monte Carlo microcanonical sampling rules. 

PACS: 75.40.Mg Numerical simulation studies 
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I — Introduction 



The practical interest of equilibrium statistical physics is to determine the canonical 
average 

<Q>T= Y. s eM-E S /T) « 

of some macroscopic quantity Q. Here, the system is kept under a fixed temperature T, 
and the Boltzmann constant is set to unity. Both sums run over all possible microstates 
available for the system, and Es (Qs) is the value of its energy (quantity Q) at the 
particular microstate S. The exponential Boltzmann factors take into account the energy 
exchanges between the system and its environment, in thermodynamic equilibrium. 
An alternative is to determine the microcanonical average 

of the same quantity Q. In this case, the energy is fixed. Accordingly, the system is 
restricted only to the g(E) degenerate microstates corresponding to the same energy level 
E, and the sum runs over them. 

The canonical average (1) can also be expressed as 

<Q>T - Z E g(E)^P(-E/T) ' (3) 

where now the sums run over all allowed energy levels E. 

Both g(E) and < Q(E) > depend only on the energy spectrum of the system. They do 
not vary for different environment conditions. For instance, by changing the temperature 
T, the canonical average < Q >t varies, but not the energy functions g(E) and < Q(E) > 
which remain the same. Canonical Monte Carlo simulations are based on equation (1), 
determining < Q >t- one needs another computer run for each new fixed value of T. 
Instead of this repetitive process, it would be better to determine g(E) and < Q(E) > once 
and forever: canonical averages can thus be calculated from (3), without re-determining 
g(E) and < Q(E) > again for each new temperature. 

The Broad Histogram method [1] (hereafter, BHM) is based on the exact relation (4), 
to be discussed later on. This equation allows one to determine g(E) from the knowledge of 
the microcanonical averages < N up (E) > and < N dn (E) > of certain macroscopic quanti- 
ties: by fixing some energy jump AE, the number N^ p (Ng n ) counts the possible changes 
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one could perform on the current microstate S, yielding an energy increment (decrement) 
of AE. Adopting some micro canonical computer simulator in order to determine these 
averages, and also < Q(E) >, one can calculate canonical averages through (3) without 
further computer efforts for different temperatures. As the Broad Histogram relation (4) 
is exact and completely general for any system [2] , the only possible source of inaccuracies 
resides on the particular microcanonical simulator chosen by the user. 

BHM does not belong to the class of reweighting methods [3-7]. These are based on 
the energy probability distribution measured from the actual distribution of visits to each 
energy level: during the computer simulation, a visit counter V(E) is updated to V(E) + 1 
each time a new microstate is sampled with energy E. At the end, the (normalized) 
histogram V(E) measures the quoted energy probability distribution. It depends on the 
particular dynamic rule adopted in order to jump from one sampled microstate to the 
next. One can, for instance, adopt a dynamic rule leading to canonical equilibrium at some 
fixed temperature To. Then, the resulting distribution can be used in order to infer the 
behaviour of the same system under another temperature T [3] . As canonical probability 
distributions are very sharply peaked around the average energy, other artificial dynamic 
rules can also be adopted in order to get broader histograms [4-7], i.e. non-canonical forms 
ofV(E). 

BHM is completely distinct from these methods, because it does not extract any 
information from V(E). The histograms for N up (E) and N dn (E) are updated to A^ up (i?)-|- 
Ng p and N dn (E) + A^| n each time a new microstate S is sampled with energy E. Thus, 
the information extracted from each sampled state S is not contained in the mere upgrade 
V(E) — > V(E) + 1, but in the macroscopic quantities A^ p and N^ n carrying a much 
more detailed description of S. In this way, numerical accuracy is much higher within 
BHM than any other method based on reweighting V(E). Moreover, the larger the system 
size, the stronger is this advantage, due to the macroscopic character of Ng P and Ng n . 

A second feature distinguishing BHM among all other methods is its flexibility con- 
cerning the particular way one uses in order to measure the fixed-i? averages < N up (E) > 
and < N dn (E) >. Any dynamic rule can be adopted in going from the current sampled 
state to the next, provided it gives the correct microcanonical averages at the end, i.e. a 
uniform visitation probability to all states belonging to the same energy level E. The rel- 
ative visitation frequency between different energy levels E and E' does not matter. Any 
transition rate definition from level E to E' can be chosen, provided it does not introduce 
any bias inside each energy level, separately. In this way, a more adequate dynamics 
can be adopted for each different system, still keeping always the full BHM formalism 
and definitions. On the other hand, multicanonical approaches [5-7] are based on a priori 



3 



unknown transition rates: they are tuned during the simulation in order to get a flat distri- 
bution of visits at the end, i.e. a uniform V(E) for which the visiting probability per state 
is inversely proportional to g(E). By measuring the actually implemented transition rates 
from E to E', which has been tuned during the computer run and must be proportional 
to g(E)/g(E'), one can finally obtain g(E). Thus, multicanonical approaches are strongly 
dependent of the particular dynamic rule one adopts. Within BHM, on the other hand, all 
possible transitions between E and E' are exactly taken into account by the quantities 
iV up and N dn themselves, not by the particular dynamic transition rates adopted during 
the computer run. 

This text is divided as follows. Section II presents the method, while some avail- 
able microcanonical simulation approaches are quoted in section III. In section IV, some 
particularities of the microcanonical ensemble are discussed. Possible improvements in 
what concerns microcanonical sampling rules are presented in section V. Conclusions are 
in section VI. 



II - The Method 



Consider a system with many degrees of freedom, denoting by S its current microstate. 
The replacement of S by another microstate S' will be denoted a movement in the space 
of states. The first concept to be taken into account is the protocol of allowed movements. 
Each movement S — > S' can be considered allowed or forbiden, only, according to some 
previously adopted protocol. Nothing to do with the probability of performing or not this 
movement within some particular dynamic process: BHM is not related to the particular 
dynamics actually implemented in order to explore the system's space of states. We need 
to consider the potential movements which could be performed, not the particular path 
actually followed in the state space, during the actual computer run. Mathematically, the 
protocol can be defined as a matrix P(S, S') whose elements are only 1 (allowed movement) 
or (forbiden). The only restriction BHM needs is that this matrix must be symmetric, 
i.e. P(S, S') = P(S', S), corresponding to microscopic reversibility: if S — > S' is an allowed 
movement according to the adopted protocol, so is the back movement S' — >■ S. Given a 
system, one needs first to define this protocol. For an Ising magnet, for instance, one can 
define that only single-spin flips will be considered. Alternatively, one can accept to flip 
any set of spins up to a certain number, or to flip clusters of neighbouring spins, or any 
other protocol. BHM does not depend on which particular protocol is adopted. This free 
choice can be used in order to improve the efficiency of the method for each particular 
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application (besides the already quoted freedom in choosing the simulational dynamics, 
which has nothing to do with this protocol of virtual movements). 

Consider now all the g(E) microstates S belonging to the energy level E [8], and 
some energy jump AE > promoting these states to another energy level E' = E + AE. 
For each S, given the previously adopted protocol, one can count the number N^ p [1] of 
allowed movements corresponding to this particular energy jump. The total number of 
allowed movements between energy levels E and E' , according to the general definition 
(2) of microcanonical average, is g(E) < N up (E) >. On the other hand, one can consider 
all the g(E') microstates S' belonging to the energy level E' and the same energy jump 
— AE < 0, now in the reverse sense. For each S' one can count the number NgV [1] of 
allowed movements decreasing its energy by AE. Due to the above quoted microscopic 
reversibility, the total number g(E') < N dn (E') > of allowed movements between energy 
levels E' and E is the same as before. Thus, one can write 

g(E) < N up (E) > = g(E + AE) < N dn (E + AE) > , (4) 

which is the fundamental BHM equation introduced in [1]. The method consists in: a) to 
measure the microcanonical averages < N up (E) >, < N dn (E) > corresponding to a fixed 
energy jump AE, and also < Q(E) > for the particular quantity Q of interest, storing the 
results in .E-histograms; b) to use (4) in order to determine the function g{E); and c) to 
determine the canonical average < Q >t from (3), for any temperature. Step a) could be 
performed by any means. Step b) depends on the previous knowledge of, say, g(0), the 
ground state degeneracy. However, this common factor would cancel in step c). 

There is an alternative formulation of BHM [9], based on a transition matrix ap- 
proach [10]. Other alternatives can be found in [11-14]. Interesting original analyses were 
presented in those references. All of them differ from each other only on the particular 
dynamics adopted in order to measure the BHM averages < N up (E) > and < N dn (E) >. 
The common feature is the BHM equation (4). For multiparametric Hamiltonians, the 
energy E can be replaced by a vector (E\, E<i . . .): in this way the whole phase diagram 
in the multidimensional space of parameters can be obtained from a single computer run 
[15], representing an enormous speed up. 

Besides the freedom of choosing the protocol of allowed movements, the user has also 
the free choice of the energy jump AE. In principle, the same g(E) could be re-determined 
again for different values of AE. Consider, for instance, the uniform Ising ferromagnet on 
a L x L square lattice with periodic boundary conditions, and only nearest-neighbour links, 
for an even L > 2. The total energy can be computed as the number of unsatisfied links 
(pairs of neighbouring spins pointing in opposed senses), i.e. E = 0, 4, 6, 8, 10 . . . 2L 2 — 6, 
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2L 2 — 4, and 21? [16]. Adopting the single-spin-nip protocol of movements, there are 
two possible values, namely AE = 2 and AE = 4, for the energy jumps. Thus, one 
can determine the same g(E) twice, within BHM. For that, one could store four distinct 
E-histograms: N up (AE = 4), N dn (AE = 4), N up (AE = 2), N dn (AE = 2). The first 
one corresponds, for each microstate, to the number of spins surrounded by four parallel 
neighbours, whereas the second to the number of spins surrounded by four neighbours 
pointing in the opposed sense: together, they determine g(E) through (4), with AE = 4, 
from the previous knowledge of g(0) = 2 and g(6) = 4L 2 . The third and fourth histograms 
correspond, for each microstate, to the number of spins surrounded by just three or just 
one parallel neighbours, respectively. Together, they can be used, with AE = 2, in order 
to determine g(E) also through (4), from the previous knowledge of g(4) = 2L 2 . 

In practice, when Monte Carlo sampling is used as the instrument to measure the 
micro canonical averages, this freedom on the choice of AE can also be used in order 
to improve the statistical accuracies [1,2,17-19], by taking all the possible values of AE 
simultaneously. Provided one has always AE « E (hereafter the ground state energy 
will be considered as E = 0), one can store only two .E-histograms, with the combinations 

N^ = J2Ws P ^E)] 1/AE (5a) 

AE 

and 

Nt = ^Wl n (m} 1/AE , (5b) 

AE 

counted at each averaging state. This trick is an approximation very useful in order to 
save both memory and time. However, it could, in principle, introduce systematic errors, 
depending on the particular application, and should be avoided in those cases. 

Ill — Some Microcanonical Simulation Approaches 

The averaged value defined by equation (2) refers to all microstates S[E] corresponding 
to the same energy level E, each one being counted just once. For large systems, their 
number g(E) is normally very large (except, in most cases, near the ground state). In 
order to obtain a Monte Carlo (random sampling) approximation for < Q(E) >, one 
actually averages only over a restricted number of microstates, much smaller than g(E). So, 
these selected averaging states must represent the whole set without any bias, besides the 
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normal statistic fluctuations. The dynamic rule must prescribe exactly the same sampling 
probability for each state. The appropriate selection of these unbiased averaging states 
within the same energy level is not an easy task. It depends on the particular system 
under study, and also on the particular set of allowed movements one adopts. There is 
no general criterion available in order to assure the uniform sampling probability among 
fixed- E microstates. 

One possible microcanonical simulation strategy is to perform successive random 
movements, always keeping constant the energy. For instance, the Q2R cellular automaton 
follows this strategy concerning Ising-like models [20]. Each movement consists in: a) to 
choose some spin at random; b) to verify whether the energy would remain the same under 
the flipping of this particular spin; and c) to perform the flip in this case. There is no 
proof that this strategy is unbiased. Numerical evidences support this possibility, although 
good averages are obtained only after very, very long transients [21]. On the other hand, 
these enormous transient times can be avoided [22] by starting the Q2R dynamics from a 
previously thermalized state, i.e. by running first some canonical steps under a well chosen 
temperature corresponding to the desired energy. One interpretation of these findings is 
the following. Combined with the single-spin-flip protocol, the fixed-energy dynamics is a 
very restrictive strategy in what concerns a fast spread over the whole set of microstates 
with energy E. Indeed, numerical evidences of non ergodicity were found [23]. Neverthe- 
less, either by waiting enormous transient times or by preparing the starting states, Q2R 
remains a possible choice for microcanonical simulator of Ising systems. However, it will 
not give good averages at all for very small energies. 

An alternative strategy is to relax a little bit the fixed-energy constraint. This idea 
was introduced [24] by allowing only small energy deviations along the path through the 
space of states: successive random movements are accepted and performed, provided they 
keep the energy always inside a small window, i.e. always within a pre-defined set of few 
adjacent energy levels. Although sampling different energy levels during the same run, the 
visits to each one are taken into account separately by storing the data in cummulative E- 
histograms. Even so, each energy level could not be completely free of influences from the 
neighbouring ones. Nevertheless, this strategy could be a good approximation provided: 
a) the energy window width AE is very small compared with the energy E itself; and b) 
the final average < Q(E) > is a smooth function of E. Again, the method does not work 
very well for very low energies, where the condition AE « E cannot be fulfilled. This 
problem could be partially minimized by adopting some smart tricks [25], although the 
very low energies would never be well described [18]. 

A third strategy is to relax completely the fixed-energy constraint, by accepting any 
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energy jump. Now, as g(E) is a fast increasing function of E, one is more likely to toss 
an energy increment than a decrement. The result of accepting any tossed movement 
would be a fast and irreversible arrival at the maximum entropy region corresponding to 
infinite temperature, sampling energies only near the maximum of g(E). In order to avoid 
this, one needs to introduce some acceptance restriction for energy-increasing movements, 
trying to get a uniform sampling along the whole energy axis. This was first introduced 
[26] within the distinct context of finding optimal solutions (energy minima) in complex 
systems. The idea is to divide all the possible movements one could perform, starting 
from the current microstate, in two classes: increasing or decreasing the energy. First, 
one of these two classes is randomly tossed. Then, a random movement belonging to 
the tossed class is performed. This strategy corresponds to an energy random walk, and 
assures a uniform sampling along the whole energy axis, on average. In spite of this very 
useful feature in what concerns the search for optimum states in complex systems, this 
strategy cannot provide correct thermal averages because different energy jumps are mixed 
together, violating the relative Boltzmann weights between different energy levels. In order 
to obtain correct thermal averages, one needs to divide the possible movements not in two, 
but in as many classes as different allowed energy jumps exist: for each fixed positive AE, 
one counts the number of increasing-energy movements (E — > E + AE) and the number 
of decreasing-energy ones (E —* E — AE), the same AE for both, in order to accept or 
not the currently tossed movement. In other words, one needs to consider precisely the 
AE-dependent BHM quantities Ng P and iV| n defined in [1]. 

Many variants of this energy random walk dynamics can be defined, the best one 
[27] being a direct consequence of the BHM equation (4) itself, as follows. In order to 
obtain a flat distribution of visits along the energy axis, one needs: a) to toss a random 
movement starting from the current microstate with energy E; b) to perform it whenever 
the energy decreases; and c) to perform it only with probability g(E) j g(E + AE) , whenever 
the energy increases (AE > being the increment). From the exact relation (4), once one 
has some previous estimate of the energy functions < N up (E) >o and < N dn (E) >o, this 
probability is also equal to 

<N dn (E + AE) >p 
P"P- <N^(E)> ' [) 

The dynamic rule proposed in [27] is then based on a two-step computer simulation. First, 
one obtains an estimate of < N up (E) > and < N dn (E) > , by any means. Then, using 
this first estimate in a further, independent computer run, one performs the above defined 
dynamics, energy-increasing movements (E — > E + AE) being accepted only according 
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to the probability (6). From this second run, one measures accurate (according to [27]) 
averages < N up (E) > and < N dn (E) > from which g(E) can be finally obtained from the 
BHM relation (4). 

This flat histogram dynamic recipe was proposed [27] in order to solve some numerical 
deviations observed in previous versions which could be viewed as approximations to it. As 
discussed in section IV, unfortunately things are not so easy. These problems are not merely 
related to particular dynamic approximate recipes, but to another characteristic of the 
system itself: the discreteness of the energy spectrum. All fundamental concepts leading 
to the microcanonical ensemble are based on the supposition that all energy changes AE 
are much smaller than the current energy E. In other words, microcanonical ensemble is 
defined by disregarding the energy spectrum discreteness. That is why conceptual problems 
appear: a) at very low energies, for any system; b) at any energy, for tiny systems. A better 
understanding of these subtle concepts cannot be obtained by simple improvements of the 
dynamic recipe: a new conceptual framework allowing to treat also discrete spectra is 
needed (and lacking). 

On the other hand, excepting for the two situations quoted in the above paragraph, 
the various approximations to the acceptance rates (6) are not so bad as supposed in 
[27], according to the evidences shown also in sections IV and V. Thus, let's quote these 
approximations which make things easier. First, instead of performing two computer runs 
in order to obtain the transition rates (6) from the first and averages from the second, 
one can perform a single one gradually accumulating < N up (E) > and < N dn (E) > in 
-E-histograms. At each step, in order to decide whether the currently tossed movement 
must be performed or not, one uses the already accumulated values themselves, by reading 
both the numerator and the denominator of (6) from the corresponding histograms at the 
proper energy channels E and E + AE. Actually, this trick was already introduced (and 
used) in the original publication [1]. This approximation will be denoted by Al. It follows 
the same lines of real-time-defined transition rates of the multicanonical sampling methods 
[5-7]. A further approximation, hereafter called A2, consists in ignoring the AE appearing 
in the numerator of (6), reading both the numerator and the denominator from the same 
energy channel E. This saves two real division operations by the current number of visits 
V(E) and V(E + AE). The third additional approximation consists in taking together 
all possible energy jumps AE through equations (5a) and (5b), hereafter called A3. This 
saves computer memory, because only a pair of histograms, one for < N up (E) > and the 
other for < N dn (E) > are needed, instead of a pair for each different possible value of AE. 
Also, this approximation saves time because the statistic fluctuations are not spread over 
many histograms, the whole available data being superimposed in only two. 
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A further yet approximation is simply to replace the numerator and the denominator of 
(6) by the current values N$ p and N^ n corresponding to the current microstate S, instead 
of reading previously averaged values. In theory, for large systems, this replacement could 
be justified if N up and iV dn are shown to be self averaging quantities, in spite of the further 
fluctuations it introduces. However, contrary to the cases Al, A2 and A3 (described in the 
last paragraph and actually tested [1,2,17-19]), this procedure does not save any computer 
time or memory, being useless in practice. 

On the other hand, the procedure [28] of counting N$ p and iVj n at the current mi- 
crostate S, without classifying them according to different values of AE, is no 
longer an approximation: it is wrong in what concerns the measurement of averages, 
because the relative Boltzmann probability for different energies would be violated. It cor- 
responds to the mistake of missing the exponent 1/AE in equations (5a) and (5b). Only 
when applied to the completely different context of finding optimal solutions in complex 
systems [26], without performing averages, this procedure is justifiable. 

The original multicanonical approach [5] can be re-formulated according to the en- 
tropic sampling dynamics [6]. The degeneracy function g(E) is gradually constructed 
during the computer run. It is based on an acceptance probability g(E)/g(E') for each 
new tossed movement from the current energy E to another value E' . After some steps, 
the whole function g(E) is updated according to the new visit trials, and so on. How- 
ever, this method is yet based exclusively on the histogram for V(E). Another possibility 
[29] is to adopt exactly this same dynamic rule in order to sample the averaging states, 
measuring also the microcanonical averages < N up (E) > and < N dn (E) > during the 
computer run. At the end, g(E) is obtained through the BHM equation (4), instead of the 
values updated during the simulation: the results [29] , of course, are much more accurate. 
Note that, in this case, exactly the same microstates are visited, i.e. the same Markovian 
chain of averaging states. In this way, the better BHM performance is explicitly shown to 
be a consequence of the more detailed description of each averaging state, compared with 
reweighting, multicanonical approaches. 

Uniform probability distributions along the energy axis may be not the best strategy. 
Being independent of the particular simulation dynamics, BHM allows one to get better 
sampling statistics in some energy regions. In [19] the Creutz dynamics [24,18] is combined 
with the energy random walk [1]. First, one chooses a window corresponding to few 
adjacent energy levels. Then, the energy random walk dynamics is applied inside this 
window up to a certain predefined number of averaging states were sampled. After that, 
the window is moved one energy level up, by including the next allowed energy level at 
right and removing the leftmost. Then, the same procedure is repeated inside this new 



10 



window up to a slightly larger number of averaging states, and so on. Reaching the critical 
region, this number is kept fixed at its maximum value, before to go down again for energies 
above the critical region. In this way, BHM allows the user to design the profile of visits 
along the energy axis, according to the numerical accuracy needed within each region. 



IV — Canonical x Microcanonical Simulations 

The equivalence between the various thermodynamic ensembles (canonical and micro- 
canonical, for instance) is a widespread belief. However, this is true only in the so-called 
thermodynamic limit iV — > oo, where iV is the number of components forming the system 
under study. For finite systems, and thus for any computer simulational approach, the 
lack of this equivalence [30] poses serious conceptual as well as practical problems. 

Canonical computer simulation approaches were very well developed since the pioneer- 
ing work [31], half a century ago. By following some precise recipes (random movements 
transforming the current microstate into the next one) a Markovian chain of averaging 
states is obtained, from which thermodynamic canonical averages are calculated. As a 
consequence of this conceptual development, particular recipes were shown to provide un- 
biased canonical equilibrium. On the other hand, microcanonical simulation has never 
attracted the attention of researchers during this same half century (with some few excep- 
tions) . That is why some fundamental concepts concerning this subject are misunderstood 
in the literature. 

The direct way to construct a fixed- E, microcanonical simulator would be to accept 
a new randomly tossed movement only if it does not change the energy. However, this 
constraint could introduce non-ergodicity problems, depending on the particular set of 
movements one adopts. For the Ising model, for instance, this problem seems to occur by 
adopting only one-spin flips [21-23]. In order to minimize it, one needs to allow more-than- 
one-spin flips [32]. However, by flipping only two spins far away from each other, after 
each whole-lattice one-spin-flip sweep, the magnetic order is destroyed below the critical 
energy [32]. Thus, the strictly fixed energy approach is problematic, and the alternative 
is to relax it, allowing some energy changes. However, this also poses troubles, once the 
equilibrium features (magnetization, correlations, etc.) of the system at some energy level 
E are not exactly the same at another level E' . By travelling from E to E' without time 
enough for equilibration, one could introduce biases from instates into E' averages: 
the multiple-energy dynamics could distort the strictly uniform probability distribution 
inside each energy level. In short, the construction of a good microcanonical simulator is 
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not a simple subject. 

If you have not enough conceptual understanding about some particular subject, a 
good idea is to resort to another similar subject for which your conceptual understanding 
is already firmly stablished. Let's define a particular microcanonical simulator by using 
well stablished canonical rules, for which the temperature T is fixed since beginning, being 
an external control parameter. The canonical average of any macroscopic quantity 
Q becomes a function of T, as < Q >t in equations (1) or (3). Let's consider some 
finite system with N components, and its average energy < E >t- Although the energy 
spectrum is discrete [8], < E > T is a continuous function of T (except for a possible 
isolated temperature where a first order transition may occur). Thus, one can tune the 
value of T in order to have < E > T coincident with some previously chosen energy level E: 
let's call this tuned temperature T(E). Then, a correct microcanonical simulator recipe 
is: a) to run some of the many known canonical recipes with fixed temperature T(E); 
b) to measure the quantity Qs of interest, for each averaging state S whose energy 
is E; c) to accumulate Qs as well as the number of visits to the particular energy level 
E, during the run; and d) to calculate the microcanonical average < Q(E) > at the end, 
by dividing the accumulated sum of Qs by the number of visits to level E. This recipe, 
of course, may not be very efficient, once one will visit many energy levels, others than 
the previously fixed value E, storing data concerning only this particular level. Also, the 
precise temperature T(E) must be known a priori, and perhaps this knowledge could be 
achieved only by performing some previous runs. Moreover, the whole process must be 
repeated for each different energy level. Nevertheless, this recipe is perfectly correct, and 
will be the basis for our reasonings hereafter. 

Table I shows the exact data for a square 4x4 lattice Ising ferromagnet, obtained by 
direct counting all the 2 16 = 65.536 possible microstates. The adopted protocol of (virtual) 
movements is the set of all possible one-spin flips. From the first two columns, the exact 
temperatures T(E) can be obtained for each energy level E, by analytically calculating 
< E >t as a function of T from equation (3), and then imposing < E >t = E. Table II 
shows the data obtained from canonical simulation with fixed T(8) = 3.02866216 (in the 
usual units where the energy corresponding to each bond is ±J, instead of or 1), 10 8 
whole lattice sweeps, and 32 independent samples. Thus, the total number of averaging 
states is 3.2 x 10 9 . For all averaged values displayed, statistical fluctuations occur at most 
on the two rightmost digits. The second column shows the actual number V(E) of visits 
to each energy level. The expected statistical relative deviations are, thus, of the order of 
V(E)~ 1 ^ 2 in each case, in agreement with the ones actually obtained (roughly 1 part in 
10 4 ). Within this statistical accuracy, the coincidence between the first line in Table II and 
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the corresponding exact values for E = 8, in Table I, is a further evidence confirming that 
the micro canonical simulator defined in the above paragraph is indeed correct: it does not 
introduce any bias besides the normal statistical numerical fluctuations. 

Tables III to V show data obtained from canonical simulations with fixed temperatures 
T(10) = 3.57199419, T(12) = 4.66862103 and T(14) = 8.33883787, respectively. These 
values were tuned in order to give the exact average energies < E >t = 10, 12 and 14, 
respectively. Note again the coincidence of the second line in Table III, the third line in 
Table IV, as well as the fourth line in Table V with the exact values presented in Table I, 
within the numerical accuracy. 

Level E = 16 is just the center of the energy spectrum, corresponding to T(16) = oo. 
In order to simulate this situation, we adopted a very high temperature, namely T = 100, 
in Table VI. Its fifth line is supposed to be compared with the corresponding exact values 
in Table I. 

The other lines in Tables II to VI also coincide with the corresponding exact values 
in Table I, within the numerical accuracy. However, these further coincidences are not 
completely trustable, because Table II, for instance, was obtained by fixing the simulational 
temperature T(8) tuned in order to give the average energy < E >t= 8. Thus, this Table 
II is out of tune for the other energy levels E ^ 8. Accordingly, Tables III, IV, V and 
VI are out of tune for energies E ^ 10, E ^ 12, E ^ 14 and E ^ 16, respectively. 
Indeed, considering for instance the simulations at fixed T(E = 14), the relative deviation 
obtained for N dn (E = 8) with AE = 2 and 4 are respectively 14 or 18 times larger than 
the expected 1 part in 10 4 . In principle, only data obtained from canonical simulations 
performed at the right temperature T(E) (i.e. the one for which < E >t = E) could 
be taken as microcanonical averages for this particular energy level E. However, the 
systematic deviations induced by taking also data corresponding to neighboring energy 
levels others than E seem to be very small. In particular, for larger systems, and near 
the critical region, justifying the simple approach of adding histograms obtained at 
different temperatures [17,9]: in those cases, the temperature is a slowly varying function 
of E (dT/dE = at the critical point, in the thermodynamic limit). 

In the thermodynamic limit, the canonical temperature can be obtained by the sta- 
tistical definition 



i = iL lim , (7) 

T de JV^oo N v ; 

where e = E/N is the energy density, and coincides perfectly with the value T(E) quoted 
before, for which < E >t = E. However, for finite systems, both the thermodynamic limit 
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N — > oo as well as the derivative limit Ae = AE/N — > cannot be performed. A paliative 
procedure is simply to forget them, transforming equation (7) in 



T m (E) 



1 



In g(E + AE) - lng(E) 
AE 



(8a) 



where the subscript m means "micro canonical temperature" defined just now for finite 
systems, and AE is the energy gap between level E and some other level above it. Of 
course, in principle, T m (E) also depends on AE, which indicates a first difference be- 
tween this and the true canonical temperature. Moreover, contrary to the real, canonical 
temperature T{E), this T m (E) is not an external parameter depending on the system's 
environment and equilibrium conditions: it is simply an alternative formulation for the 
energy spectrum g(E) itself, being thus environment-independent. Note also that T m (E) 
is defined only at the allowed energies belonging to the discrete spectrum, whereas the 
real temperature T(E) (the one for which < E >t = E) can be defined for any value E, 
continuously. An alternative formula is 



The most famous canonical recipe [31] is: a) to toss some random movement, starting 
from the current state; b) to calculate the energy variation AE this movement would 
promote if actually implemented; c) to perform the movement, whenever AE < 0, counting 
one more step; and d) to perform the movement only with the Boltzmann acceptance 
probability exp(— AE/T), if AE > 0, counting one step anyway. Normally, one takes 
a new averaging state after N successive steps (one MC step). Let's stress that, here, 
the temperature T is the real, canonical one: in order to use this recipe to measure 
micro canonical averages at some fixed energy E, one must take T = T{E), the temperature 
for which < E >t= E. Instead, by taking T = T m {E) in equation (8a), it is an easy 
exercise to show that the Boltzmann acceptance probability exp(— AE/T) would be equal 
to g(E)/g(E + AE). But this is just the acceptance probability adopted within the flat 
histogram dynamics [27]. 

Figure 1 shows T m (E) calculated from equation (8a), for a 4 x 4 square Ising ferro- 
magnet. The exact values of g(E) were used. The open circles correspond to AE = 2, and 
the diamonds to AE = 4. The continuous line is the exact T(E), also calculated from the 
exact values of g(E). The deviations between T m (E) and T(E) are very strong. Note that 
there is no approximation at all, neither in T m (E) nor in T(E). The deviations represent 
true differences between canonical and micro canonical ensembles, which are indeed very 



T m (E + AE/2) 



1 



lng(E + AE) -lng(E) 
AE 



(8b) 
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strong for this tiny system. Obviously, the condition AE « E is not fulfilled, and the 
discreteness is inevitable along the whole energy spectrum. 

Figure 2 reports the same data as figure 1, with the same symbols, now for a 32 x 32 
lattice, and using equation (8b). The exact values for g(E) were taken from [33]. The 
same strong deviations occur again, but now only at the very beginning of the energy 
spectrum, where the condition AE « E does not hold, as can be seen in the upper inset. 
However, near the critical region, the deviations become much smaller, as can be seen in 
the lower inset where the vertical scale is 100 times finer than the upper one. For E > 64, 
the energy spectrum discreteness can be neglected within a very good approximation, even 
with N = 1024 being still very far from the thermodynamic limit in equation (7). The 
larger the system size, the better becomes the situation, because the relation AE « E 
becomes more and more fulfilled. However, at the very beginning of the spectrum this 
relation will never be fulfilled, even for large systems. 

Both limits in equation (7), namely the thermodynamic one N — > oo and AE — > 
corresponding to the energy derivative, were neglected in equation (8a). However, only the 
latter seems to have disastrous consequences when one uses T m (E) as an approximation 
for T(E). Indeed, even for very small systems like the 32 x 32 square lattice (not so 
tiny as 4 x 4), the deviations are very small, provided the condition AE « E holds. 
In other words, the deviations between the microcanonical temperature T m (E) and the 
true canonical value T(E) comes almost exclusively from the breakdown of the condition 
AE « E, not from finite size effects. Even in the thermodynamic limit, T m (E) 
will differ from T(E) near the ground state, due to the discretness of the 
energy spectrum. On the other hand, to take T m (E) instead of T{E) is a very, very 
good approximation far from the ground state: we have seen before that microcanonical 
averages obtained from canonical simulations are not sensitive to temperatures out of tune 
for a given energy, moreover when the deviations are of the order of that shown in the 
lower inset of figure 2. 

Tables II to VI present good numerical accuracies not only for the tuned tempera- 
tures, but also when out-of-tune values for T were considered. Moreover, based on the 
reasonings above, this feature will be even improved for larger and larger systems. The 
imediate consequence is the possibility of adding different histograms for < N np (E) > and 
< N dn (E) >, obtained from distinct canonical simulations with different temperatures, as 
already tested in [17,9]. Of course, this approach is much more efficient than to take a 
different canonical simulation with fixed temperature T(E) for each different energy level 
E, without superimposing the histograms. However, it still needs many computer runs, 
one for each fixed temperature. 
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In order to improve even more the efficiency, one can try the following strategy. First, 
one defines the Boltzmann acceptance probability exp[—AE/T(E)] for each energy level 
E and each possible energy jump AE. Note that this is not the same as canonical simu- 
lations where the acceptance probability exp(— AE/T) depends only on AE but not on 
E. Then, by following this non-canonical, .E-dependent acceptance probability, one runs a 
single computer simulation, visiting the whole energy axis, accumulating the histograms 
for < N up (E) > and < N dn (E) >. This approach is similar to the flat histogram dy- 
namics [27], using the true temperature T(E) instead of the microcanonical value T m (E). 
Table VII shows the results for the same tiny system already considered before. Surpris- 
ingly, strong deviations appear. In order to analyse the source of these deviations, let's 
introduce another very similar alternative, adopting a different acceptance probability 
exp[— AE/T(E + AE/2)]. The canonical Boltzmann probability is taken at the center 
E + AE/2 of the interval correponding to the energy jump from E to E + AE, not at 
the current energy E. This symmetrization trick is supposed to diminish the numerical 
deviations due to impossibility of performing the limit AE — > 0, i.e. the derivative in 
equation (7), for this tiny system. The results are shown in Table VIII, in complete agree- 
ment with the exact results, Table I. Thus, the source of the deviations in Table VII is 
the lack of the condition AE « E. By taking into account both the current energy level 
E as well as the (would-be) next energy E + AE in a symmetric way, the numerical 
deviations were eliminated. Of course, for larger systems and far from the ground state, 
where the condition AE « E holds, the differences between Tables VII and VIII would 
also disappear. 

The flat histogram dynamics uses another acceptance probability, namely 

where E and E + AE also play a symmetric role. Simulational results are shown in Table 
IX, where the exact values for g(E) (or, alternatively, < N up (E) > and < N dn (E) >) were 
adopted in order to determine the acceptance probabilities (9). The results are again 
coincident with the expected ones, Table I [34]. On the other hand, if one tries to break 
the quoted symmetry, ignoring AE in the numerator of the right-hand-side term in (9), 
i.e. approximation A2, numerical deviations similar to Table VII would appear, for this 
tiny system. 

As one does not know a priori g(E) (or, alternatively, < N up (E) > and < N dn (E) >), 
one can use some previous estimates < N up (E) > and < N dn (E) > , and adopts the 
acceptance probability (6). In order to follow this recipe [27], one needs to determine the 
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quoted estimates from a previous computer run. According to [27], these previous esti- 
mates do not need to be very accurate. However, actual numerical tests [29] show results 
for < N up (E) > and < N dn (E) > worse than the inputs < N up (E) > and < N dn (E) > 
themselves! A better possibility is the random walk dynamics originally used in order to 
test the broad histogram method [1]. It is the same as the flat histogram, equation (6) or 
(9), with the approximation Al (optionally, also A2 and A3) quoted in last section. Al 
consists in taking the current, already accumulated values of the histograms N up (E) and 
N dn (E + AE), in real time during the computer run itself, instead of the true averages at 
the right-hand side of equation (6) or (9). Results for the same tiny system treated before 
are shown in Table X. This approach is the same as the multicanonical real-time-measured 
transition probability already adopted in other earlier methods, for instance the entropic 
sampling [6] . A criticism to this approach is its non-strictly-Markovian, history-dependent 
character. According to [29], it is actually better than the fixed transition probability pro- 
posed in [27]. Approximation A2 consists in neglecting AE in the numerators of equation 
(6) or (9), breaking the symmetry between levels E and E + AE, and cannot be applied 
to this tiny system due to the spectrum discretness. This approach A2 is also shown to 
violate a particular detailed balance condition [27]. A3 consists in taking all possible val- 
ues for AE together, by using equations (5a) and (5b). Both approximations A2 and A3 
(but not Al) are bad: a) at very low energies, for any system; b) at any energy, for tiny 
systems. 

The difference between the dynamics adopted in Tables VII to X and the true canonical 
rule adopted in Tables II to VI is the following. Canonical simulations adopt the same 
acceptance probability exp(— AE/T) for any tossed movement increasing the energy by 
AE, no matter which is the current energy E. As a consequence, the energies visited during 
the run became restricted to a narrow window around the canonical average < E >t- 
The larger the system size, the narrower is this window. On the other hand, within the 
non-canonical rules adopted in Table VII to X, the acceptance probabilities depend also 
on E (and E + AE). As a consequence, instead of a narrow distribution of visits, one 
gets a broad distribution covering the whole energy axis. This feature is, of course, a 
big advantage over canonical rules, because only one run would be enough to cover a 
large temperature range. Nevertheless, the numerical results could be wrong, depending 
on the actual dynamic rule one adopts (Table VII, for instance). A better conceptual, 
theoretical understanding of these and other independent dynamic rules is needed, and 
lacking. Concerning BHM, the only constraint to be considered is the uniform probability 
visitation inside each energy level, separately. For other, reweighting methods based on 
the actual visitation profile V(E), also the detailed relative distribution between different 
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energy levels must be considered. 

Up to now, we have tested many different dynamic rules in order to measure the 
micro canonical averages < N up (E) > and < N^ n (E) > from which one can determine the 
desired quantity g(E) by BHM equation (4). The most efficient approaches are the in- 
dependent rules (Tables VIII, IX and X), where a single computer run is enough. Among 
them, under a practical point of view, the random walk dynamics [1] corresponding to 
Table X is the best choice, once one does not need any previous knowledge about the 
quantities < N up (E) > and < Nd n (E) > to be measured. This approach corresponds to 
approximation Al. The other two further approximations A2 and A3 [1] improve even 
more the efficiency. However, due to the energy spectrum discreteness, they are limited 
by the constraint AE « E: one needs to avoid them for tiny systems, or very near the 
ground state even for large systems, where this constraint cannot be fulfilled. All this 
matter corresponds to the subject discussed in reference [27]. Let's now discuss another, 
subtle, further possible source of unaccuracies, which may appear when one abandon the 
safe canonical dynamical rules and adopts non-canonical, independent dynamics in order 
to sample the whole energy axis during the same computer run. 

In order to introduce the subject, let's resort again to canonical simulations, where 
some temperature value T is fixed since beginning. Imagine one starts such a simulational 
process from a randomly chosen microstate S : its energy E$ would be in general far from 
the average value < E >t, as well as many other features of this microstate which would 
be far from their equilibrium counterparts. By plotting the energy of each successive 
averaging microstate as a function of the time, one would get a curve fluctuating around 
an exponential decay to the final value < E >t- at the end, only statistical fluctuations 
around this constant value remain. Thus, the very beginning of the Markovian chain of 
states will give wrong (biased) contributions to the final averages < Q >t- These are the 
so-called thermalization transient steps. Nevertheless, they represent no problem at all, 
because one can always take an enormous number of microstates after this transient bias 
is already over, pushing the systematic deviations to below any predetermined tolerance. 
Better yet, one can simply discard the contribution of these initial states, by starting the 
averaging procedure only after the transient is over. In the statistical physicists' jargon, one 
can assert that "the system is already thermalized" , after these initial out-of-equilibrium, 
transient steps. 

As quoted in [27], detailed balance is a delicate matter. Detailed balance conditions 
are useful only in order to ensure that the final distribution, chosen by the user, will be 
the correct one, for instance the Gibbs distribution for equilibrium canonical averages. 
However, these conditions do not ensure this final distribution will be reached within a 
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finite time. In our case of interest, i.e. the independent, broad-energy, non-canonical 
dynamics, the system never reaches canonical equilibrium inside each energy level. All 
these dynamic rules correspond to generalized Boltzmann factors exp[— AE/Q(E)], where 
the "temperature" 0(E) varies from one microstate to the next, along the Markovian 
chain. Concerning the energy, for instance, instead of an exponential decay to the canon- 
ical equilibrium situation, one gets a random walk visiting all energies, with fluctuations 
covering the whole energy axis all the time. The same large fluctuation behaviour holds 
also for other quantities, in particular the one for which the microcanonical averaging is in 
progress. This dynamics follows an eternal transient in what concerns the safe canonical 
framework. This feature may introduce systematic numerical deviations. Although those 
possible out-of-canonical-equilibrium problems cannot be observed in our Tables VIII, IX 
and X, they could be crucial for larger systems where large energy jumps in few steps be- 
come possible: features of a particular far energy region which the system recently cames 
from could introduce biases in the current energy averages. In this eternal-transient case, 
detailed balance conditions and all the related theorems give little help. To construct a 
good, efficient microcanonical simulator, by visiting the whole energy axis during a single 
computer run, seems to be a much more delicate matter. How to assure a uniform prob- 
ability distribution inside each energy level also covering broad regions of the energy 
spectrum? This problem is open to new insight, new ideas. 

Exemplifying how difficult would be to analyse the uniformity of visits within each 
energy level, let's take a simple example: a L x L square lattice Ising magnet (L > 4, 
N = L 2 spins), with E = 8. Considering the magnetization density m, level E = 8 is 
divided into three classes. The first one contains gi(8) = N(N — 5) states with only two 
non-neighbouring spins pointing down (all other spins up, or vice- versa), corresponding to 
m = 1 — 4/N. The second class consists of (72(8) = 12 N states with a 3-site cluster of 
neighbouring spins pointing down, with m = 1 — 6/N. The third class has (73(8) = 2N 
states presenting a plaquette of four spins pointing down, and m = 1 — 8/N. No single- 
spin flip can transform a state with E = 4 or E = 6 into another state belonging to this 
third class, for instance. In order to assure a uniform distribution of visits, this lack must 
be compensated by other possible single-spin flips from E = 10 and E = 12. On the 
other hand, the visitation frequency to each one of these higher-energy states depends on 
feeding rates from higher yet energies, and so on. One must prescribe adequate transition 
probabilities to each such movement, taking into account all its consequences on the next, 
next-next step, and so on. Chessboard is an easier game. Fortunately, the first class 
containing ~ A 2 states dominates the counting for large lattice sizes. For L = 32, for 
instance, gi(8) = 1,043,456 states are in the first class, representing 99% of the whole 
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number g(8) = <7i(8) + (72(8) + (73(8) = 1, 057, 792. Similar behaviours also occur for higher 
energies, giving us a solace: we remain with the hope that single-spin flips may lead to the 
desired uniformity for large enough systems. For a tiny 4x4 lattice, on the other hand, 
things go worse, once two further classes must be added to E = 8 level: (74(8) = 16 states 
with a single line of spins down, corresponding to m = 1/2; and (75(8) = 8 states with 
two neighbouring lines of spins down, with m = 0. Then, only (71(8) = 176 states out of 
(7(8) = 424, i.e. 42%, belong to the first class. This feature, of course, partially explains 
the bad results in Table VII. However, a detailed explanation for the good results obtained 
in Tables VIII, IX and X, following the same reasonings, is not easy. 



V - Improved Microcanonical Simulators 

While we have not yet a perfect and efficient microcanonical simulator, the advan- 
tages of the Broad Histogram method, i.e. the exact and completely general equation 
(4), are already in hands. The only problem is to obtain good estimates for the averages 
< N up (E) > and < N dn (E) >, as functions of E. One paliative solution is to use the avail- 
able, not-yet perfect microcanonical simulators, perhaps with some repairing procedures. 
Even without any repair, the random walk naive approach introduced in [1] can be useful. 
First, observing Tables II to VI, one notes that some (large) degree of inaccuracy on the 
temperature can be tolerated still within very good precision on the final results. Second, 
comparing figures 1 and 2, one can note that larger systems present a plateau of almost 
constant temperatures covering a large energy range. This range is just the interesting 
one, namely the critical region where the specific heat diverges for N — > 00. Third, with 
small variations of the temperature from one microstate to the next, within this region, 
the out-of-canonical-equilibrium problems should be strongly diminished. 

Anyway, we can try to introduce some repairing procedures into the random walk 
dynamics [1] (or into any other independent transition rate dynamics) . The first idea is to 
equilibrate the system before measuring any averaging quantity at the current microstate. 
This can be acomplished by running some canonical steps just before the measuring proce- 
dure. Suppose one gets some microstate with energy E, during the random walk dynamics. 
As this microstate comes from other energy levels, submitted to acceptance probabilities 
others than the currently correct value exp[— AE/T(E)], it is supposed to carry some un- 
desired biases. Then, one can simply include some canonical steps with fixed temperature 
T = T{E), in order to let the system relax to .E-equilibrium, before taking the averages. 
T(E) can be estimated at each step from the current, already accumulated histograms for 
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N up (E) and N dn (E) (approximation Al), and the further approximations A2 and A3 can 
also be adopted. All these tricks were already introduced in the original publication [1] . 

Another repairing procedure is simply to forbid large energy jumps. The simplest 
way to perform this task is to count a new averaging state after each single-spin-flip trial. 
Normally, one adopts N trials, i.e. a whole lattice sweep, before counting a new averaging 
state, in order to avoid possible correlations along the Markovian chain of states. In our 
micro canonical case, however, different energy levels correspond to independent averaging 
processes, and one can try to abandon this precaution. Of course, this approach also saves 
a lot of computer time. Table XI shows the results obtained for a 32 x 32 lattice, where 
approximation Al (real-time measurement of the transition rates) was adopted, with 10 9 
averaging states per energy level. The expected relative error due to finite statistics is thus 
3 x 10 -5 . Indeed, the observed deviations coincide with this, in spite of the figure 0.008 
[27] predicted by detailed balance arguments. Even adopting the further approximations 
A2 (which explicitly violates detailed balance) and A3, the errors remain the same (last 
column). Once more, the possible source of unaccuracies has nothing to do with detailed 
balance dictating the relative frequency of visits between different energy levels. On 
the contrary, the only crucial point is the uniformity of visits inside each energy level, 
separately. 

Perhaps the procedure of measuring averages after each new single-spin-flip trial could 
indeed introduce some bias, through some undesirable correlations along the Markovian 
chain, although this possibility is not apparent at all in Table XL Perhaps this bias could 
appear only for much larger statistics. An alternative way to avoid large energy jumps still 
taking averages only after each whole lattice sweep is to restrict the random walk to narrow 
energy windows. This corresponds to the Creutz energy-bag simulator [24] combined with 
the random walk dynamics [1] inside each window. It was introduced in [19], and the 
results also show relative errors much smaller than those predicted by detailed balance 
arguments. The specific heat for a 32 x 32 square lattice is shown in figure 3, where 
approximations Al, A2 and A3 were used. By restricting the energy to narrow windows, 
one is also restricting the temperature to small variations, thus forcing the system to be 
always near the canonical equilibrium conditions for any energy inside the current window. 
Once again, the best performance occurs at the critical region. 



VI - Conclusions 

The Broad Histogram Method (BHM) introduced in [1] allows one to determine the 
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energy spectrum of any system, i.e. the degeneracy g(E) as a function of the energy E, 
through the exact and completely general equation (4). First, one needs to adopt some 
reversible protocol of allowed movements in the system's space of states. Reversible means 
that for any allowed movement S — > S' , the back movement S' — > S is also allowed. For 
Ising models, for instance, one can choose single-spin flips, to flip clusters of neighboring 
spins, etc. In fact, one could invent any protocol. For each state S, N^ p counts the 
number of such allowed movements for which the system's energy would be increased by a 
fixed amount AE. Accordingly, Ng n counts the number of allowed movements decreasing 
the energy by the same amount AE. The energy jump AE is also chosen and fixed since 
the beginning. < N up (E) > and < N dn (E) > are the microcanonical, fixed- E averages 
for these two quantities, both functions of the energy. Then, in order to determine g(E) 
from the BHM equation (4), one needs only to find a way, any way, to measure those 
microcanonical averages. 

By following the same way adopted in measuring < N up (E) > and < N dn (E) >, one 
can also obtain the microcanonical average < Q(E) > of the particular thermodynamic 
quantity Q of interest (magnetization, density, correlations, etc). All these microcanonical 
averages are independent of the particular environment the system is currently interact- 
ing with. In other words, they do not depend on temperatures, equilibrium conditions, or 
any other thermodynamic concept: they are determined by the system's energy spectrum 
alone. Thus, after g(E) is already determined through the BHM equation (4), once and 
forever, the same system can be submitted to different environment conditions, and its 
behaviour (equilibrium or not) can be studied resorting always to the same g(E). Within 
the particular case of canonical equilibrium, for instance, the thermal average < Q >t of 
the quantity Q can be obtained from equation (3), for any temperature T. If one decides 
to use computer simulations as the instrument measuring < N up (E) >, < N dn (E) > and 
< Q(E) >, then only one computer run is enough to determine the whole temperature 
dependence, continuously, without need of repeating again the simulation for each new 
temperature. 

Other computer simulation methods, the so-called multicanonical sampling strategies 
[5-7], also allow the direct determination of g(E). All of them are based on the counting 
V(E) of visits to each energy level E. Every time each energy level E is visited by the 
current state S, along the Markovian chain, the counting is updated from the current V(E) 
to V(E) + 1. Within BHM, however, the ^-histograms for A^ up and A^ dn are updated from 
the current N up (E) and N dn (E) to N up (E) + N^ p and N dn (E) + Nj n , respectively, for 
the same current state S. This corresponds to macroscopic upgrades instead of the mere 
counting of one more state. Thus, a much more detailed information is extracted from each 
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averaging state withing BHM than any other method, giving rise to much more accurate 
results. Moreover, the larger the system size, the stronger is this advantage, due to the 
macroscopic character of the BHM quantities N up and N dn . 

Another advantage of BHM over other methods is its complete independence concern- 
ing the particular dynamic rule adopted in order to measure the micro canonical averages. 
Any recipe can be used, provided the correct i?-functions < N up (E) >, < N dn (E) > 
and < Q(E) > were accurately determined. As the fundamental equation (4) is exact, 
any numerical deviation observed on the final results are due to the particular measuring 
recipe one chooses to adopt, not to the method itself. This is a big advantage, once one 
can choose a more adequate measuring instrument for each different system under study. 
The only requirement to obtain correct micro canonical averages is a uniform probability 
distribution inside each energy level, separately. Detailed balance between different 
energy levels is irrelevant for BHM. All possible transitions between different energy 
levels are exactly counted by the BHM quantities N up and N dn themselves, not by the 
particular stochastic dynamic recipe one adopts. On the other hand, besides the uniform 
sampling probability inside each energy level, multicanonical methods also depend on the 
particular transition rates between different energies, which are tuned during the computer 
run in order to get a flat distribution of visits at the end. Detailed balance is fundamental 
for multicanonical methods. Thus, besides the acccuracy advantage commented before, 
any good dynamic rule for multicanonical sampling is also good for BHM, but the reverse 
is not true. Indeed, BHM allows the user to design his own profile of visits along the 
energy axis, taking for instance a better statistics near the critical region. 

During the past half century, theoretical studies provided us with many recipes con- 
cerning detailed balance for Markovian processes, leading to Gibbs equilibrium distribu- 
tion. These recipes show us how to design adequate transition rates between different 
energy levels. Unfortunately, there are no equivalent theoretical studies in what concerns 
adequate rules leading to a uniform probability inside each energy level. Perhaps this miss- 
ing point is due to the fact that most studies were directed towards canonical distributions, 
where a sharp region of the energy axis is visited: after some transient steps (normally 
discarded from the Markovian chain) the system becomes trapped into this sharp region. 
Thus, the spread over this restricted region is supposed to occur naturally, giving rise 
to the required uniformity inside each level. At least for smooth independent quantities 
< Q(E) >, the quoted sharpness solves the problem. However, this is not the case if one 
needs to cover wide energy ranges. For both multicanonical sampling as well as BHM, this 
is an open problem waiting for new ideas, new insight. How to assure a uniform probability 
of visits inside each energy level? Some clues towards the answer were discussed, and some 
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very efficient practical solutions were proposed in section V. 
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Tables 



E 


9(E) 


< N up (E) > 


< N up (E) > 


< N t[e (E) > 


< N dn (E) > 


< N dn (E) : 






(AE = 4) 


(AE = 2) 


(AE = 0) 


(AE = 2) 


(AE = 4) 





2 


16 














4 


32 


11 


4 








1 


6 


64 


8 


6 





2 





8 


424 


6.11321 


6.33962 


1.81132 


0.90566 


0.83019 


10 


1728 


3.55556 


7.03704 


3.55556 


1.55556 


0.29630 


12 


6688 


2.31101 


5.97129 


5.51196 


1.81818 


0.38756 


14 


13568 


1.13208 


5.58491 


5.88679 


2.94340 


0.45283 


16 


20524 


0.75307 


3.69207 


7.10973 


3.69207 


0.75307 



Table I Exact data concerning the 4x4 square lattice Ising ferromagnet. The second half of 
the spectrum (E = 18 . . . 32) is symmetric the to first half (E = 14 ... 0). N tie is the 
number of potential movements (one-spin flips) keeping the energy unchanged. 



E 


V(E) 


< N up (E) > 


< N up (E) > 


< N tie (E) > 


< N dn (E) > 


< N dn (E) > 




xlO 8 


(AE = 4) 


(AE = 2) 


(AE = 0) 


(AE = 2) 


(AE = 4) 


8 


5.62 


6.1145 


6.3376 


1.81155 


0.90606 


0.83028 


10 


6.11 


3.5565 


7.0366 


3.5540 


1.55630 


0.29665 


12 


6.32 


2.3111 


5.9711 


5.5119 


1.81843 


0.38744 


14 


3.42 


1.13151 


5.5851 


5.8879 


2.9429 


0.45259 


16 


1.38 


0.75302 


3.6916 


7.1109 


3.6913 


0.75315 



Table II Canonical simulation for the 4x4 square lattice Ising ferromagnet, with a fixed 
temperature T(8) = 3.02866216 for which the exact energy average is < E > T = 8. 
For all averaged values, statistical fluctuations fall on the two last digits, at most. 
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E 


V(E) 


. a Tim / t — 1\ ~. 

< N up (E) > 


. a Tim / 7 — ?\ ^ 

< N up (E) > 


< N tie (E) > 


< N an (E) > 


< N an (E) > 




xlO 8 


/ A 7~~ 1 A \ 

(AE = 4) 


(AE = 2) 


(AE — 0) 


( A 7~~ 1 C\ \ 

(AE = 2) 


/ A 7~~ 1 A \ 

(AE = 4) 


8 


4.68 


6.1123 


6.3410 


1.81107 


0.90577 


0.82989 


10 


6.22 


3.5559 


7.0367 


3.5552 


1.55578 


0.29639 


12 


7.85 


2.3113 


5.9713 


5.5115 


1.81836 


0.38771 


14 


5.20 


1.13148 


5.5858 


5.8870 


2.9426 


0.45307 


16 


2 57 


0.75274 


3.6914 


7.1116 


3.6917 


0.75257 


! Ill 


The same 


as Table II, now with a fixed temperature T(10) = 3.57199419 for which 




the exact 


energy average 


is < E > T = 


10. 






E 1 


V(E) 


< N up (E) > 


< N up (E) > 


< N t[e (E) > 


< N dn (E) > 


< N dn (E) > 




xlO 8 


(AE = 4) 


(AE = 2) 


(AE = 0) 


(AE = 2) 


(AE = 4) 


8 


3.08 


6.1122 


6.3404 


1.8116 


0.90667 


0.82909 


10 


5.34 


3.5550 


7.0372 


3.5552 


1.55529 


0.29597 


12 


8.78 


2.31065 


5.9707 


5.5138 


1.81777 


0.38711 


14 


7.56 


1.13206 


5.5846 


5.8872 


2.9436 


0.45255 


16 


4.85 


0.75310 


3.6923 


7.1091 


3.6924 


0.75303 



Table IV The same again, now with a fixed temperature T(12) = 4.66862103 for which the 
exact energy average is<i?>T=12. 



E 


V(E) 


< N up (E) > 


< N up (E) > 


< N t[e (E) > 


< N dn (E) > 


< N dn (E) > 




xlO 8 


(AE = 4) 


(AE = 2) 


(AE = 0) 


(AE = 2) 


(AE = 4) 


8 


1.30 


6.1137 


6.3404 


1.8098 


0.90444 


0.83169 


10 


3.27 


3.5550 


7.0379 


3.5554 


1.55540 


0.29628 


12 


7.84 


2.31063 


5.9716 


5.5123 


1.81816 


0.38735 


14 


9.83 


1.13243 


5.5845 


5.8866 


2.9436 


0.45291 


16 


9.20 


0.75357 


3.6916 


7.1094 


3.6920 


0.75339 



Table V The same once more, now with a fixed temperature T(14) = 8.33883787 for which the 
exact energy average is < E >t = 14. 
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8 
10 
12 
14 
16 

Table VI 



E 

8 
10 
12 
14 
16 

Table VII 



E 

8 
10 
12 
14 
16 

Table VIII 



V(E) 


. a Tim / t — 1\ 

< N up (E) > 


. a Tim / 7 — 1\ ^ 

< N up (E) > 


< iV tie (.E) > 


< N an (E) > 


< N an (E) > 


xlO 8 


/ A 7~~ 1 A \ 

(AE = 4) 


(AE = 2) 


(AE = 0) 


( A 7~~ 1 C\ \ 

(AE = 2) 


(AE = 4) 


0.35 


6.1132 


6.3405 


1.8106 


0.90464 


0.83108 


1.40 


3.5560 


7.0366 


3.5552 


1.55577 


0.29642 


5.22 


2.31103 


5.9717 


5.5112 


1.81838 


0.38769 


10.17 


1.13244 


5.5847 


5.8863 


2.9436 


0.45298 


14.77 


0.75265 


3.6922 


7.1102 


3.6922 


0.75266 



Again, now with a fixed temperature T = 100, mimicking T(16) = oo. 



V(E) 


< N up (E) > 


< N up (E) > 


< N t[e (E) > 


< N dn (E) > 


< N dn (E) > 


xlO 8 


(AE = 4) 


(AE = 2) 


(AE = 0) 


(AE = 2) 


(AE = 4) 


3.91 


6.1075 


6.3168 


1.8413 


0.93716 


0.79731 


4.11 


3.5647 


7.0105 


3.5752 


1.55947 


0.29019 


4.96 


2.32037 


5.9336 


5.5420 


1.83372 


0.37031 


3.99 


1.13602 


5.5671 


5.8935 


2.9675 


0.43584 


3.55 


0.73739 


3.6917 


7.1232 


3.7290 


0.71873 



Random walk, pseudo-canonical simulation with asymmetric -B-dependent Boltzmann 
acceptance probability exp[— AE /T(E)], to be compared with the exact values in 
Table I. 



V(E) 


< N up (E) > 


< N up (E) > 


< N t[e (E) > 


< N dn (E) > 


< N dn (E) > 


xlO 8 


(AE = 4) 


(AE = 2) 


(AE = 0) 


(AE = 2) 


(AE = 4) 


3.25 


6.1134 


6.3392 


1.81132 


0.90608 


0.82998 


3.87 


3.5571 


7.0343 


3.5561 


1.55637 


0.29609 


5.48 


2.31315 


5.9677 


5.5123 


1.81968 


0.38715 


5.57 


1.13256 


5.5833 


5.8871 


2.9458 


0.45130 


6.56 


0.75220 


3.6918 


7.1109 


3.6939 


0.75118 



Random walk, pseudo-canonical simulation with symmetric independent Boltzmann 
acceptance probability exp[— AE /T(E+AE /2)], to be compared with the exact values 
in Table I. 



29 



E 


V(E) 


. a Tim / t — 1\ ^ 

< N up (E) > 


. a Tim / 7 — ?\ ^ 

< N up (E) > 


< N tie (E) > 


< N an (E) > 


< N an (E) > 




xlO 8 


/ A 7~~ 1 A \ 

(AE = 4) 


(AE = 2) 


/ A 7 — ' r\ \ 

(AE = 0) 


(AE = 2) 


/ A 7~~ 1 /i \ 

(AE = 4) 


8 


4.53 


6.1123 


6.3407 


1.81151 


0.90579 


0.82973 


10 


4.53 


3.5553 


7.0376 


3.5552 


1.55561 


0.29630 


12 


4.53 


2.3112 


5.9704 


5.5130 


1.81807 


0.38733 


14 


4.53 


1.13190 


5.5848 


5.8873 


2.9432 


0.45272 


16 


4.53 


0.75283 


3.6920 


7.1108 


3.6912 


0.75321 



Table IX Flat histogram [27] simulation with symmetric E-dependent Boltzmann acceptance 
probability exp[— AE /T m (E)], equation (9), to be compared with the exact values in 
Table I. 



E 


V(E) 


< N up (E) > 


< N up (E) > 


< N tie (E) > 


< N dn (E) > 


< N dn (E) > 




xlO 8 


(AE = 4) 


(AE = 2) 


(AE = 0) 


(AE = 2) 


(AE = 4) 


8 


3.51 


6.1131 


6.3399 


1.81121 


0.90558 


0.83023 


10 


4.65 


3.5553 


7.0372 


3.5561 


1.55532 


0.29619 


12 


17.03 


2.3106 


5.9717 


5.5124 


1.81788 


0.38749 


14 


3.81 


1.13212 


5.5850 


5.8865 


2.9434 


0.45296 


16 


2.13 


0.75298 


3.6917 


7.1108 


3.6915 


0.75305 



Table X Random walk dynamics [1] simulation with approximation Al (see text), to be com- 
pared with the exact values in Table I. 
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E 



onn 

300 


1. 


.749921 


1 


.749897 


1 


.749800 


1 


.749992 


302 


1 


.748650 


1 


.748599 


1 


T AO A 0£? 

.748496 


1 


T A oa TO 

.748678 




1 


.747398 


1 


.747387 


1 


.747288 


1 


.747365 


60b 


1 


.746164 


1 


.746151 


1 


.746055 


1 


i-r a r* -i a a 

.746144 




1 


rr a ac\ An 

.744946 


1 


.744960 


1 


.744834 


1 


>~7 A A (\(\ r 

.744905 


310 


1. 


.743743 


1 


.743748 


1 


.743638 


1 


.743660 




1 


.742554 


1 


7,1 nr nn 

.742590 


1 


.74z45z 


1 


7/10/100 

.742438 


O "1 A 

314 


1 


1"7 A ~\ n/7r 

.741375 


1 


.741424 


1 


.741305 


1 


.741192 


316 


1 


.740206 


1 


.740295 


1. 


.740148 


1 


.740000 


318 


1. 


.739045 


1. 


.739128 


1 


.738995 


1. 


.738762 


320 


1 


.737889 


1. 


.737970 


1. 


.737828 


1. 


.737596 


322 


1. 


.736737 


1 


.736831 


1. 


.736655 


1. 


.736406 


324 


1 


.735586 


1. 


.735680 


1 


.735494 


1. 


.735275 


326 


1 


.734434 


1. 


.734560 


1. 


.734342 


1 


.734086 


328 


1 


.733279 


1 


.733420 


1. 


.733174 


1 


.732904 



Table XI Data for the 32 x 32 square lattice, near the critical region. Exact values are in 
the second column, and the next 2 columns correspond to different BHM runs. The 
random walk dynamics [1] with approximation Al (see text) was adopted. A new 
averaging state is counted after each single-spin flip, instead of waiting for a whole 
lattice sweep. Nearly 10 9 states per energy were sampled for ech run. The same 
numerical accuracy was obtained also by using the further approximations A2 and A3, 
last column, in spite of the detailed balance violation [27] between different energies, 
forced by A2. Within BHM, one does not need to care about the relative balance 
of visitations to different energies. Only a uniform sampling probability inside each 
energy level, separately, is important. 



31 



Figure Captions 

Figure 1 Exact canonical temperature as a function of the energy (continuous line), for a 4 x 4 
square lattice Ising ferromagnet. Symbols correspond to the microcanonical version of 
the temperature, equation (8a), with energy gaps AE = 2 (circles) or 4 (diamonds). 

Figure 2 The same as figure 1, now for a 32 x 32 lattice, with the same symbols, and equation 
(8b). Note that the deviations between the true canonical temperatures and the mi- 
crocanonical versions are now restricted to the very beginning of the energy spectrum 
(upper inset). At the critical region these deviations are much smaller (lower inset, 
with a 100 times finer scale). The deviations become smaller yet for larger systems. 

Figure 3 Specific heat for the 32 x 32 lattice Ising ferromagnet, adopting the random walk 
dynamics within restricted energy windows [19]. Circles show the exact values [33]. 
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